program diffraction_and_energy;
uses common;
uses tools,inp,outp,ppot,decor,d04_header,diffract,relax;

Var
   i,ii,j,k,l,m,n,o,kk,ll    : integer;
   iters,nsteps,mcs,nte	     : integer;
   txt,txpars,txrep	     : text;
   tinp,dmcs,cmcs,hloops     : integer;
   cmp,epot,rf,rwf,chin,chix : myreal;
   rfint,rwfint		     : myreal;
   ueff,rate,rmsd	     : myreal;
   se			     : rsites;
   dse			     : rsites3;
   w			     : v3;
   txout		     : text;
   domc			     : boolean;
   ss			     : mystring;
   chi3,echi3		     : v3;
begin
   readln(tinp);
   readln(ticoname);
   export_occup:=false;
   reset(txt,'INPUT/objects.def'); {need to set repre5, but also dpmodel}
   readln(txt,deconame);
   readln(txt,o); {nu_offset}
   readln(txt,dpmodel);
   readln(txt,o); if(o>0) then repre5:=true else repre5:=false;
   close(txt);
   BASIS;
   reset(txt,ticoname);
   IMPORT_TILING_BODY(tinp,txt);
   TABLE_INIT;
   if (tinp>0) then begin
      IMPORT_CONFIG(txt);
   end;
   close(txt);
   UNIT_CELL(multi);  
   {for i:=1 to 3 do begin
   for o:=1 to 3 do blt[o,i]:=bl3[i,o];
   end;}
   PP_INP;
   if (tinp=0) then begin
      i:=0;
      READ_DECO(i); {reads objects.def, ap.def}
      for i:=1 to tnodes do nan[i]:=0;
      for i:=1 to tnodes do NEW_NODE_OBJECTS(i,1);
   end;
   nat:=0; for i:=1 to tnodes do
      for j:=1 to nan[i] do if (n_a[i,j,1]>0) then nat:=nat+1;
   if debug[13] then begin
      writeln('Atomic decoration: ',nat:1,' atoms.');
   end;
   {note: ch_status needed for relax, but its called in dfit_init}
   nntable(rcplus);
   UPPAV(epot,se);
   if debug[13] then writeln('E_ini = ',epot/nat:1:8,' eV/at');
   debug[131]:=true;
   READ_DFIT_INP(dmcs,mcs,hloops,nte);
   DFIT_INIT;
   domc:=((mcs>0) or (dmcs>0));
   save_e:=1;
   if(dmcs>0) then rewrite(recu,'udf.rep');
   rewrite(txrep,'dfit.rep');
   writeln(txrep,'#TL HL Uef E RWF RF chi_dat chi_dw chi_e');
   rewrite(txpars,'dfit.copas');
   export_tiling_header(txpars);
   export_tiling_body(txpars);
   DFIT_START_INFO(txpars);
   for kk:=1 to nte do begin
      for ii:=1 to hloops do begin
	 nntable(rcplus);
	 if (domc) then begin
	    if(dmcs>0) then begin
	       DCMC(dmcs,kk,epot,ueff,rf,rwf,chin);
	       etarg:=ueff;
	    end;
	    if(mcs>0) then begin
	       CMC(mcs,kk,cueff,sigmadr,rate,epot,ueff,rmsd);
	       etarg:=ueff;
	    end;
	 end else begin {converting beta to energy}
	    lambdafit:=lambdas[kk];
	    etarg:=etargs[kk];
	 end;
	 for o:=1 to 3 do echi3[o]:=0;
	 for ll:=1 to nfita do begin
	    if ((fita[ll,1]<0) or (fita[ll,2]<0)) then begin
	       if (fita[ll,1]<0) then randomize_params(1); {positions}
	       if (fita[ll,2]<0) then randomize_params(0); {dw factors}
	    end else begin
	       lambdafit:=lambdas[kk];
	       etarg:=etargs[kk];
	       DFIT(fita[ll,1],fita[ll,2],ueff,rf,rwf,chi3);
	       if(fita[ll,1]>1) then echi3:=chi3;
	       chin:=chi3[1];
	    end;
	 end;
	 chi3[3]:=echi3[3];
	 chin:=chi3[1];
	 if(cueff>1e-10) then UPPAV(epot,se) else epot:=ueff;
	 paradump(txpars,kk,ii,epot,ueff,rf,rwf,chi3);
	 if(chin<chimin)then begin	    
	    chimin:=chin;
	    HKL_OUT('DFIT/hkl.min');
	    ss:='Exporter: FitMin ';
	    DFIT_OUT('DFIT/copa.min',ss,epot,ueff,rwf,rf,chi3);
	    DPARS_OUT(ss,epot,ueff,chin,rwf,rf);
	 end;
	 writeln(kk:1,' ',ii:1,' ',ueff/natoms:1:5,' ',epot/natoms:1:5,' ',chin:1:3,' ',rwf:1:4,' ',rf:1:4,' #Ueff E CH2 RWF RF');
	 write(txrep,kk:1,' ',ii:1,' ',ueff/natoms:1:5,' ',epot/natoms:1:5,' ',rwf:1:4,' ',rf:1:4);
	 for o:=1 to 3 do write(txrep,' ',chi3[o]:1:3); writeln(txrep);
      end;
   end;
   {final call to update}
   INTERNAL_RWFAC(rf,rwf,chix);
   writeln('Approximant-internal wR R: ',rwf:1:4,' ',rf:1:4);
   HKL_OUT('hkl.out');
   ss:='Exporter: EndFit ';
   DFIT_OUT('copa.dfit',ss,epot,ueff,rwf,rf,chi3);
   if(dmcs>0) then close(recu);
   writeln(txpars,'end_paradump');
   export_tiling_atoms(txpars);
   close(txpars);
   close(txrep);
end.
